A two-species model of a two-dimensional sandpile 
surface: a case of asymptotic roughening 
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Abstract We present and analyze a model of an evolv- 
ing sandpile surface in (2 + 1) dimensions where the dy- 
namics of mobile grains (p(x, t)) and immobile clusters 
(/i(x, t)) are coupled. Our coupling models the situation 
where the sandpile is flat on average, so that there is no 
bias due to gravity. We find anomalous scaling: the ex- 
pected logarithmic smoothing at short length and time 
scales gives way to roughening in the asymptotic limit, 
where novel and non-trivial exponents are found. 

PACS 05.10.-a, 64.70.qj, 68.35.Ct 



1 Introduction 

The study of sandpile surfaces unifies the field of granu- 
lar physics [1] with that of surface growth [2] , which lat- 
ter has been the subject of prolonged theoretical analy- 
sis. One of the best known equations in surface growth 
is the Edwards- Wilkinson equation [3 , which has been 
used across a swathe of fields: from the experimental 
analysis of the epitaxial growth of thin films [4] to the 
theoretical study of anomalous growth exponents in the 
presence of drift [5 ,6 . 

However, its use in the study of sandpile surfaces 
has been limited by its one- variable formulation, viz. 
the response of surface height fo(x, t) to the applica- 
tion of stochastic perturbations. Computer simulations 
of shaken granular packings in three dimensions |7] and 
cellular automaton models [8 suggested that a full de- 
scription of the dynamics of granular surfaces needed 

Bandan Chakrabortty 

S N Bose National Centre For Basic Sciences 

Block- JD, Sector-Ill, Salt Lake, Kolkata 700098, India 

E-mail: bandan@bose.res.in 

Anita Mehta 

S N Bose National Centre For Basic Sciences 

Block- JD, Sector-Ill, Salt Lake, Kolkata 700098, India 

E-mail: anita@bose.res.in 



the coupling of granular clusters Q and flowing (mobile) 
grains. Since the growth and decay of clusters in a given 
region directly affects the surface height, it was thought 
appropriate [9] to couple the surface height /i(x, £), with 
a variable p(x, t) denoting the density of flowing grains 
across the surface. Equations were accordingly written 
down [9lHQ] which embroidered the Edwards- Wilkinson 
equation with a variety of (usually) nonlinear trans- 
fer terms coupling /i(x, t) with p(x, £), reflecting dif- 
ferent physical situations such as pouring/shaking a 
flat/sloping sandpile. The response to such perturba- 
tions could involve, for example, grains leaving their 
clusters and beginning to flow (/i(x, t) decreasing and 
p(x, t) increasing at (x, t)) \ or a current of flowing grains 
getting lodged in a dip on the surface (p(x, t) decreas- 
ing and /i(x, t) increasing at (x, t)). The purpose of the 
transfer term was to model such exchanges between the 
two species, and in so doing ensure the conservation of 
grains. A lively discussion ensued about the nature of 
these transfer terms, and many different ones were pro- 
posed to match evolving surfaces under different phys- 
ical conditions [TTj [T2 , 13 , 14 , 15 . Subsequently, the ap- 
plication of such noisy nonlinear two-species equations 
was extended to even more diverse situations, ranging 
from avalanching to ripple and dune formation p ^ fTT l 

he]. 

All of the above approaches were unified by one par- 
ticular feature: sandpile surfaces exhibited simple scal- 
ing, independent of the length and time scales consid- 
ered. However, experiments suggested that anomalous 
scaling might arise [19] . with different exponents char- 
acterising short and long lengths/timescales. In par- 
ticular, it was suggested that sloping sandpile surfaces 
might exhibit asymptotic smoothing under tilt [20] ; sim- 
ulations suggested that this might be due to the large 
avalanches released, which 'swept the surface clean' [2T] . 

1 We use clusters to indicate collections of grains on sur- 
faces which are relatively immobile 
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New coupled equations (again drawing on the Edwards- 
Wilkinson equation [3]) were accordingly written down 
in one dimension, which explicitly investigated the ef- 
fect of tilting the surface of a sloping sandpile: they 
showed that although roughening continued to be ob- 
served at short length and time scales, sandpiles mani- 
fested asymptotic smoothing [22 . The presence of bias 
is crucial to such smoothing, of course - avalanches 
flow down, rather than up a sloping sandpile. It was 
then natural to ask: what happens in the absence of 
such a bias, i.e. when a sandpile is flat on average? 
Also, what happens in dimensions higher than one? The 
present paper asks, and then answers, these questions: 
anomalous scaling is observed once again, but this time 
around, the crossover is from logarithmic smoothing for 
short length/time scales to asymptotic roughening. 

We now review some pertinent facts. In general, re- 
sults of even single-species equations in more than one 
dimension, are hard to obtain analytically. An excep- 
tion is the Edwards- Wilkinson equation [3], which can 
be solved in any space dimensions (d): it is critical for 
d = 2, with roughening exponents (to be defined be- 
low) a = 2j3 = 1 — d/2, where logarithmic scaling is 
observed [2]l23]. Even in the case of the KPZ equation 
[24] where an analytical solution exists for d = 1 [25j 
[26] , one has to resort to approximate numerical solu- 
tions and compare them with conjectures in d = 2 [27] 

We end this introductory section with a recap of 
scaling laws pertinent to surface growth [30,31 . In gen- 
eral, power-law scaling is manifested by growing sur- 
faces, so that the interfacial width W(L, t) of an initially 
flat interface continues to grow as a power of the time 
t until it saturates: at saturation, its width W sa t(L,t) 
grows with a power of the system size L. Thus: 



W 2 (L,t) ~t 
W? at (L,t)~L 2 



2/3 



(1) 

(2) 



where f3 is the growth exponent and a is the rough- 
ness exponent. These two asymptotic power laws can 
be combined into a single scaling form, 



W\L,t)~L 2a f(t/L z ), 



(3) 



with z = a/P the dynamic exponent and / an appro- 
priate scaling function. 

However, power-law scaling is sometimes replaced 
by logarithmic scaling, especially for surfaces at their 
critical dimensionality [2 . In such cases, we have: 



W 2 (L,t) 



' Ht/a), 
In(L /a), 



(4) 
(5) 



where a is a lattice constant. These two asymptotic ex- 
pressions can be combined into a single scaling form 
[231132] to give: 



W 2 (L,t) - In 



\ L t 




- 9 1 
a y 





(6) 



where g is an appropriate scaling function. 

The plan of this paper is as follows: In Section [2] we 
present and describe our model equations; our results 
are presented and analysed in Section [3} We discuss our 
results and make some concluding remarks in Section 

Efl 



2 The model equations 

Our model equations in (2 + 1) dimensions contain a 
coupling between /i(x, t), the local height of immobile 
clusters and p(x, t), the density of mobile particles; this 
is done through a transfer term T, as in previous ap- 
proaches [1 , and we follow closely the notation used 
earlier [TP] : 

dh(x,t) 



dt 
dp(x, t) 
dt 



= V^(x,£)-T + 7^(x,£), 
= V 2 p(x,t)+T + ^(x,£), 



(7) 
(8) 



-p(x, t) V 2 /i(x, t) + (u- /xp(x, t))\ V/i(x, t) | 



The transfer term models exchanges between immo- 
bile clusters and flowing grains: the first part is propor- 
tional to surface diffusion, while the second part de- 
pends on the magnitude of the local slope at any point. 
This term is thus symmetric with respect to the sign of 
the slope, and suggests that any deviation from a flat 
slope results in the conversion of grains from one to the 
other species. In more detail: 

The first term p(x, t)V 2 /i(x, t), models the contribu- 
tion to surface diffusion made by the 'effective' surface 
of flowing grains when there is continuous avalanching 
(p(x, t) 7^ for all x and t) across the sandpile. The sec- 
ond term z/|Vft(x, t)\ represents the spontaneous gener- 
ation of flowing grains in response to (any direction of) 
tilt. The term /ip(x, t)|V/i(x, t)\ represents the effect 
of a boundary layer below which the surface dynam- 
ics will not penetrate by limiting the release of flowing 
grains to this cutoff depth. The introduction of this 
term can be seen to provide a regulator, so that tilt- 
ing will only ever be a moderate perturbation. After a 
short while when the condition (dp(x,i)/dt) = ob- 
tains, then the average value of p(x, t) will be given by 
(p(x, t)) ^ v I \i. This condition, apart from its physical 
reasonableness, is essential for numerical simulations, 
as it ensures that /i(x, t), p(x, t) and their fluctuations 
remain finite and measurable. Thus fi plays the role of 
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an experimental/numerical cut off [331 134] and the ratio 
of v I [i controls the inherent dynamics. 

The effects of shaking a sandpile (randomising the 
cluster distributions /i(x, £)) and pouring grains on a 
sandpile (randomising the distribution of p(x,t)) are 
modelled by the noise terms r/^(x, t) and r/ p (x, t) re- 
spectively. These are accordingly chosen to be two in- 
dependent Gaussian white noises, characterised by their 
widths Ah, A p as follows: 

( % (x, t)r, h &, t')) = A 2 h 6 d (x - x')6(t - f ), (9) 
(t? p (x, t) Vp ( X ', f )> = A 2 p 6 d (x - x')6(t - t'). (10) 

3 Numerical results and analysis 

The simulations have been performed by discretizing 
Equation [7] and Equation [8| both in space and time. 
For this purpose we have implemented multigrid algo- 
rithms [35] on a two-dimensional square lattice. Since 
it is critical behaviour that interests us, we have chosen 
relatively large spatial grids, (a — Ax = Ay ~ 1.0); in 
order to avoid numerical instabilities, we have chosen 
relatively small time stpng, At ~ 0.001. We have used 
a regulator function F(T) to eliminate numerical diver- 
gences from the transfer term T as in earlier work [lOj 
l22l. This ensures that: 



-1 for T < -1 
T for -1 < T < 1 
1 for T > 1 



F(T) = 

We start with the following intuition: 



(ii) 



— Given the comments made in the introduction, we 
would expect the unperturbed behaviour of h(x,t) 
to be logarithmically smooth, since the Edwards- 
Wilkinson equations has d c = 2. 

— We would likewise expect critical fluctuations in p(x, t) 
to be manifest only for v ^> /i, because the mean 
value of p(x, t) is constrained by the relation (p(x, t)) ~ 
vjp (see above). 

Before evaluating critical exponents, we examine the 
real-space behaviour of both species as a function of the 
relevant parameters. According to the analytical pre- 
dictions [2 ,23] mentioned above, we expect W 2 (L,t) to 
grow logarithmically with time for d = 2, i.e W 2 (L, t) ~ 
hit /a: accordingly, in Figure [l] we plot W 2 (L,t) vs 
hit /a for different values of /i, v. Although we have 
checked that our results do not depend on system size 
by checking for systems of size (32 x 32, 64 x 64, 128 x 
128,256 x 256), we present, for conciseness, only our 
results for L = 64 in Fig. [I] To begin with, in Figure 
[TJl, we take v = 0; our findings confirm the intuition 
that for any (non-zero) value of /i, the /i(x, t) species 





Fig. 1 Square widths W 2 (L, t) are plotted vs hit /a with 
L = 64 for both h and p species and for different val- 
ues of /i, v. Black indicates the h species and the red 
indicates the p species in every graph, (a) A represen- 
tative plot for p ^ (chosen to be p = 1.0, in this 
case) and v = 0. The qualitative behavior of the h — p 
curves does not depend on the precise value of p with 
in a large range, (b) The plot for p = v = 0. (c) When 
both p and v are non-zero, the fluctuations of h and 
p become comparable at a special point (see text). For 
p = 1.0, this occures at v s = 4.7; the corresponding 
graphs are presented, (d) The system size dependence 
of the coupling parameters at the special point (/i s , v s ) 
is explored by plotting (p Sl v s ) vs. 1/L. 



is critical, while p(x, t) has fully non-critical fluctua- 
tions. As soon as we set p = 0, keeping v = 0, we 
find - surprisingly - that the critical fluctuations of 
p(x, t) begin to dominate those of /i(x, t) (Figure [TJd) . 
In order to identify the point at which the critical fluc- 
tuations of both species become comparable, we fix 
p = 1.0 and tune v\ the results for this special point 
are shown in Figure [l]j^] We examine the system size 
dependence of this special point by analysing systems 
of size 32 x 32, 64 x 64, 128 x 128, 256 x 256; the criti- 
cal behaviour of both species is unchanged if we rescale 
the values of p and v with system size. This rescaling is 
made explicit in Figure [l]i, which suggests these special 

2 The initial bump in these plots for h(x,t) and p(x, t) re- 
flects the large initial value of p(x, t) used in order to ensure 
the positivity of p(x,f) throughout; this however is not ex- 
pected to affect the long-time scaling behaviour, as our plots 
confirm 
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values ji s and v s depend systematically on system size 
as ~ 1/L. 

In order to evaluate the critical exponents, we eval- 
uate the structure factors S(t sat , k) and S(L,uS), which 
are related to the interfacial width W 2 (L,t) as follows 



w 



\L,t) = ^S{n,t) = Ws{L,^). 



(12) 



L d ^ V,J t d 

K CO 

Because of the above relation, the scaling relations Equa- 
tion [l] and Equation [2] can be recast as, 

K~ 2 ~ 2a («->()), (13) 
S(L,u) ~uj- 2 ~ 2 P (a; -►()), (14) 
here, t sat is the saturation time of the interfacial width 
W(L,t) and the wave vector k 



We mention first of all that the critical behavior of 
the Edwards- Wilkinson equation in 2d is logarithmic 
(a = P = 0). In the work of this paper, crossovers 
to different nontrivial exponents are manifested as a 
function of different parameters. In particular, we em- 
phasizes that all the nontrivial scaling behaviors found 
here are anomalous, i.e. shows different scaling for long 
length and time scales. 

We proceed by first evaluating the cases considered 
in Figure [T] corresponding to the physical situation 
where there is noise in both species (i.e. A\ = A 2 > 0). 
The plots of structure factor for each set of parameter 
values considered are presented successively, while the 
exponent values are compiled and tabulated in Table 



la and Table lb (see below). 

The first case is presented in Figure [2iJ for the pa- 
rameter values v = and \i ^ 0. Based on Figure [TJi, 
we would not expect to see nontrivial critical fluctua- 
tions for p(x, £), an expectation borne out by Table [Ta| 
and Table [lb] Physically, this is because the absence of 
tilt [y = 0) halts the spontaneous generation of mobile 
grains in the steady state, when < p(x, t) >~ 0; the 
only source of flowing grains is ?7 p (x, t). 

The case of v — \i = is examined in Figure |2n} 
where Figure [TJd would lead us to expect that the fluc- 
tuations of p(x, t) would dominate those of /i(x, t); a 
perusal of Table [Ta| and Table [lb] confirms that this is 
the case. In this case, the dynamical equation for ft(x, t) 
can be viewed as a linear Edwards- Wilkinson equation 
having an effective diffusion constant of (1 + p(x,£)). 
The special point (/i s ,^ s ) characterises the crossover 
point where the fluctuations of p(x, t) just begin to be- 
come 'dangerous', i.e. critical, and of the same order as 




In co 



(i) Coupling parameters used were: v = 0.0, /i = 4.0. 




(ii) Coupling parameters used were: v = 0.0, /i = 0.0. 




(iii) Coupling parameters used were: v s = 4.7, /i s — 1.0. 




(iv) Coupling parameters used were: v = 15.0, fi = 1.0. 

Fig. 2 Plots of structure factors for the different cases 
(in the presence of noise in both the species: A\ = A 2 = 
1.0): (a) hi S (t sa t, k) vs In ft plot and (b) \nS(L,uj) vs 
lno; plot. 



Above the special point (fi s fixed, v > v 3 ), we would 
expect that the fluctuations of p(x, t) would dominate 



those of h(x, t). The plots of Figure [2m] which lead to those of h(x, t) (as in the case of Figure [2h]); our results, 



the exponents tabulated in Table [Ta] and Table [Tb] con- 
firm this intuition: both species show nearly comparable 
asymptotic roughening. 



presented in Figure [2rv] and tabulated in Table [Ta] and 
Table |lb[ confirm this. This can be understood intu- 
itively by realising that for the p(x,t) equation (and 
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when p(x,t) is finite), the term p(x, i)V 2 h cannot be 
renormalised away as an addition to an effective diffu- 
sion constant, as it can in the case of the /i(x, t) equa- 
tion. 

Finally, we study the situation when there is only 
sustained noise in p(x, t): i.e. A\ = after a transient 
and A 2 > 0. From a study of the model equations in one 
dimension under these conditions [10 , we would expect 
that the fo(x, t) landscape would freeze up after a tran- 
sient (with possibly non-trivial roughening exponents), 
while the p(x, t) landscape would diffuse over this frozen 
landscape with Edwards- Wilkinson exponents. We find 
indeed (Figure [3j Table [Ta| and Table lb ) that after 
a transient, the fo(x, t) landscape is frozen, with non- 
trivial roughening exponents entirely inherited from the 
transient noise in fo(x, t). At this point, the p(x, t) land- 
scape gets decoupled from the fo(x, t) landscape, and 
shows logarithmic roughening. 




Fig. 3 Plots of structure factors in the absence of noise 
in h species (A\ = after a transient and A 2 > 0): (a) 
In S(t sat , k) vs In ft plot and (b) \nS(L,uj) vs \iiuj plot. 
The coupling parameters used were: v = 15.0, /i = 1.0 
but the results hold for arbitrary values. 



Before leaving this section, we summarise our find- 
ings: except when the tilt term vanishes, the critical 
fluctuations of p(x, t) dominate those of fo(x, t) as long 
as there is noise in both species A\ > 0, A 2 > 0. 
When, however, A\ = after a transient and A 2 > 0, 
the p(x, t) species is logarithmically smooth, while the 
frozen fo(x, t) landscape is characterised by a novel non- 
trivial roughening exponent, entirely inherited from the 
transient phase. 



4 Discussion 

The main focus of our paper was the investigation of 
a two-species model of sandpile dynamics, representing 
the height profile and the dynamics of granular flow 
along a fiat surface in the presence of perturbations 
such as shaking and/or pouring. Importantly, no bias 



Mode of coupling 


Species 


a 


for large k, 


for small k 


Edwards- Wilkinson 


/i(x, t) 


0.0 


0.0 


z/ = 0,/i^0 
A\ = A* p >0 


/i(x, t) 
p(x,t) 


0.00 ±0.05 
-0.02 ±0.04 


0.62 ±0.06 
-0.02 ±0.04 


v = 0,p = 
Al = A*>0 


h(x,t) 
p(x,t) 


0.07 ±0.08 
0.00 ±0.05 


0.07 ±0.08 
0.81 ±0.07 


fi s (special), u s > p s 
Al = Aj>0 


h(~K, t) 
p(x,t) 


0.00 ±0.03 
0.00 ±0.06 


0.66 ±0.06 
0.85 ±0.05 


v,n(> v s ,ii s ),v > n 
Al = Al>0 


/i(x, t) 
p(x,t) 


0.00 ±0.05 
0.00 ±0.04 


1.28 ±0.04 
1.85 ±0.06 


arbitrary p, v 
Al = 0,Al>0 


h(x,t) 
p(x,t) 


1.45 ±0.08 
0.00 ±0.04 


1.45 ±0.08 
0.00 ±0.04 


(a) Values of a. 


Mode of coupling 


Species 


P 


for large uj 


for small uj 


Edwards- Wilkinson 


/i(x, t) 


0.0 


0.0 


v = 0,p ^ 
Al =A^>0 


/i(x, t) 
p(x,t) 


0.00 ±0.03 
-0.03 ±0.04 


0.50 ±0.04 
-0.03 ±0.04 


u = 0,fi = 
Al =A*>0 


h(x,t) 
p(x,t) 


0.01 ±0.05 
0.00 ±0.04 


0.01 ±0.05 
0.76 ± 0.06 


u s , lis (special), u s > /i s 
Al = Al>0 


h(~K, t) 
p(x,t) 


0.00 ±0.04 
0.00 ±0.06 


0.64 ±0.02 
0.76 ±0.04 


V, /i(> V s , /JL S ) , V > /Jl 

Al = Al>0 


h(x,t) 
p(x,t) 


0.00 ±0.05 
0.00 ±0.08 


0.72 ± 0.05 
0.85 ±0.07 


arbitrary /i, v 
Al = 0, A 2 p > 


h(x,t) 
p(x,t) 


0.16 ±0.03 
0.00 ±0.07 


0.16 ±0.03 
0.00 ±0.07 



(b) Values of (3. 

Table 1 Numerically estimated values of a and /3, for 
the various physical situations described in this paper. 



exists in this system, so that avalanches of grains can 
flow in both directions: this is to be contrasted with 
the downward flow of grains due to the biasing field 
of gravity, on a sloping sandpile. Given that the latter 
scenario had led to asymptotic smoothing [21,22 , we 
wanted to investigate the effect of making our exchange 
terms fully symmetric, to model a fiat sandpile surface. 

The particular model that we have explored has an 
exchange term that includes the effect of tilt, as well as 
that of a finite boundary layer, below which the bulk 
of the sandpile is protected from the dynamics at the 
surface. A major difference from the investigation of 
similar equations in one dimension [10] is that anoma- 
lous scaling is obtained in general for the most rele- 
vant species: wherever appropriate, short- wavelength 
logarithmic smoothing behaviour at small wavelengths 
crosses over to asymptotic roughening. This is largely 
due to the action of the tilt term, which in this unbiased 
case, causes flowing grains to be generated independent 
of the sign of the local slope: dips and bumps alike ac- 
cumulate flowing grains, a change from the asymmetric 
(biased) situation considered in [22] . 
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We now discuss our results in detail. When v = 0, 
flowing grains are not generated from clusters; the only 
source for these is from the 'pouring' source, rj p (x.,t). 
The form of the transfer term then implies that there 
is a net transfer from the flowing grains to the clus- 
ters, resulting in an enhanced roughness as the cou- 
plings become relevant; the /i(x, t) species is then char- 
acterised by a crossover from logarithmic smoothing to 
roughening with the novel exponents presented in the 
second line of each of Tables [Ta| and [lb] , while the 
p(x, t) species remains logarithmically smooth as ex- 
pected. Such a simplification also occurs (sixth line of 



Tables la and lb) when after a transient period when 



both noises are present, only A 2 remains finite. After a 
transient period, the /i(x, t) profile is frozen into a state 
of extreme roughness, fed by the flowing grains (whose 
dynamics show the expected EW logarithmic smooth- 
ing): the large roughening exponent obtained is novel, 
especially because it is independent of parameter val- 
ues. Finally, for v > fi such that {y y ' \i) > (z/ s //i s ), the 
critical fluctuations of p(x, t) dominate those of /i(x, t), 
which can be viewed as a state of continuous avalanch- 
ing, when the effective surface is defined by the film 
of flowing grains [22 . Anomalous scaling is again ev- 
erywhere observed, with a crossover from logarithmic 
smoothing to asymptotic roughening: this can be given 
an interpretation in the context of 'thick avalanches' 
[T5] , so that the novel roughening exponents we have re- 
ported for p(x, t) can be analysed in greater depth with 
respect to earlier studies [14 . They can also be mea- 
sured using sophisticated visualisation techniques on 
traditional rotating cylinder experiments, which have 
been the subject of considerable recent activity [36] . 
and which can fruitfully be employed in studies of the 
kind of surface roughening reported here. 
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